Nearest-neighbour Attraction Stabilizes Staggered Currents in the 2D Hubbard Model 
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Using a strong-coupling approach, we show that staggered current vorticity does not obtain in the 
repulsive 2D Hubbard model for large on-site Coulomb interactions, as in the case of the copper oxide 
superconductors. This trend also persists even when nearest-neighbour repulsions are present. However, 
staggered flux ordering emerges when attractive nearest-neighbour Coulomb interactions are included. 
Such ordering opens a gap along the (7r,0) — (0, tt) direction and persists over a reasonable range of 
doping. 



Doped Mott insulators such as the cuprates are rid- 
dled with competing ordering tendencies. In addition to 
antiferromagnetism, dx2-y2 pairing, and charge ordering, 
staggered orbital antiferromagnetism produced by local 
circulating currents has joined the list of physical states 
that might occupy a central place in the phase diagram of 
the cuprates. Namely, staggered flux order ^ might be 
at the heart of the elusive psuedogap phase. It has been 
proposed that such a state produces a genuine gap in the 
single-particle spectrum at the (tt, 0) point and competes 
with d-wave pairing leading to the termination of super- 
conductivity in the underdoped regime. Physically, cur- 
rents which alternate in sign from plaquette to plaquette 
in the copper-oxygen plane comprise the staggered flux 
phase. Consequently, the staggered flux phase can be 
thought of as a density wave having d^2_y2 symmetry. 
As charge density wave and superconducting order are 
well known to compete in strongly-correlated systems, a 
psuedogap arising from a d-density wave state has imme- 
diate appeal. 

The earliest theoretical work [|[ which showed that 
staggered flux states might be the relevant low-energy 
excitations of 2D Heisenberg antiferromagnets was based 
on 1/iV expansions of the t — J model. Such expan- 
sions cannot access the strong-coupling limit relevant to 
the copper oxides. Numerical simulations j^] using a 
variational Gutzwiller-projected d-wave superconducting 
state for the t — J model on a 10 x 10 lattice revealed 
power law decay of the staggered current. Leung Q in- 
troduced two holes into a 32 site t-J cluster and also 
found nonzero staggered vorticity correlations. More re- 
cently, density matrix renormalization group (DMRG) 
calculations j5| on two-leg ladders have found the rung- 
rung current correlation function to decay exponentially 
indicating an absence of staggered flux ordering. Conse- 
quently, convincing computational evidence for staggered 
flux ordering is lacking. In addition, there has been no 
study of the staggered flux phase based on a model in 
which double occupancy is explicitly retained. 

To fill this void, we start with the Hubbard model. 
We use the Hubbard operator approach as these opera- 



tors ^,0] are tailor-made for treating the strong-coupling 
regime, U ^ t because they diagonalize the on-site 
Coulomb repulsion. We have shown recently Q that this 
approach yields a stable d^2_y2 superconducting state 
in the absence of dynamical self-energy corrections. We 
show here that an equivalent treatment of the d-density 
wave state leads to an absence of any long-range order 
of this type when U is at least on the order of the band- 
width (U — 8t), the relevant magnitude for the cuprates. 
Absence of d-density wave ordering at this level of theory 
is significant because the dynamical or quantum correc- 
tions tend to destroy long-range order. Hence, we con- 
clude that for U — 8t, d-density wave ordering is absent 
at any filling from the 2D Hubbard model. However, 
we find that when a nearest-neighbour Coulomb interac- 
tion, V = —0.25t, is present, the staggered flux phase is 
stabilized. A gap opens along the (tt, 0) — (0, tt) direc- 
tion as anticipated in Ref. (1). Our results suggest that 
Coulomb attractions, perhaps arising from phonons 
could be relevant to the phenomenology of the cuprates. 

The starting point for our analysis is the Hubbard 
model 



U ^ nil nil -I- ^ VijUiUj 



(1) 



with tij — taij and Vij = Vaij where aij = 1 when 
ij are nearest-neighbours and zero otherwise and U the 
on-site Coulomb repulsion. For the cuprates, U ~ 
and t — Q.heV. As U is the largest energy scale in 
the problem, we use the eigen operators of the atomic 
limit, namely r^i^ = CicrUi^a and = Ci„{\ - Ui^a), 
which deflne the upper and lower Hubbard bands, respec- 
tively. These operators diagonalise the on-site interaction 
term, Uni-^nii = ^ Via''li(T hence are the natural 
starting point for a strong-coupling analysis. The Hub- 
bard operators are cumbersome in that they do not obey 
usual fcrmionic anticommutation relationships. However, 
equation of motion techniques have proven quite useful 
in circumventing the lack of a diagrammatic scheme . 
We provide here only a sketch of the method as the full 
details have been published elsewhere 0. First, we break 
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the system into two sublattices and hence consider the 4- 
component spinor, 



defined on two sites with 

\ li(T / \ li'a 



(2) 



(3) 



The key quantity of interest is the corresponding Green 
function, 



G(z, ^^ t; j, f, t') = e{t - t'){{i:„{i, i' , t),i^^ {j, f , t'}) 



(4) 



where {X, Y} is the anticommutator and (• • •) is the ther- 
mal average. Thermal averages involving operators on 
two sites 



(5) 



will in general contain a real part, {AB') and in the flux 
phase, an imaginary part, {AB')i as a consequence of 
the local broken time-reversal symmetry. Here, 7ij = 1 
if ij are nearest neighbours along the x-axis, whereas 
7y = — 1 if the two neighbours are along the y-axis and 
zero otherwise, Q = (tt, tt) so that the imaginary term 
alternates in sign between the two sublattices and the 
prime superscript signifies that the operators A and B 
are defined on nearest neighbour sites. 

To develop a system of self-consistent equations, we 
start by writing the equations of motion, 



,■(2) 



for the Hubbard operators in the presence of a finite 
chemical potential, fi. The new operator tt is defined 
by 



The remaining components, j'^^^ and j*^^' are constructed 
simply by interchanging A and B and i and i' in j^^^ and 
respectively. In the static (or pole) approximation 
§, the Gr een function equations are closed using the 
Roth B projector 



E 

In 



(8) 



which explicitly removes the components of j*^*) which are 
orthogonal to the Hubbard basis. The overlap matrix I = 



F.T.{{iljcr{i,i'),'il}l{j,j')}) is an enitirely diagonal 4x4 
matrix with 111 = 133 = 1 — (n_cr) and /22 = ^44 = {n-a). 
Here F.T. signifies the Fourier transform. At this level of 
theory, the energy levels are sharp as the Green function 



G(k,w) = (tjl-E)-i/: 
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(9) 



has poles at the real energies, Ei{k), where i — 1,2 
indexes the lower Hubbard band and i = 3,4 the up- 
per Hubbard band. Because of the two sublattice con- 
struction, i?i(k) = £^2(k + Q) and likewise, i?3(k) = 
E4k + Q). InEq. {§, a, = z{E,) / ]\^^^{E^- Ej) where 
z{uo) = Det(wl - E)(u;l - The poles of the Green 

function result from diagonalising the matrix 

E = Ml^ (10) 

where andM = F.T.{{ja{i,i'),ipl{i,j')}) isa4x4matrix 
defined on the two sublattices. The matrix elements of M 
contain two types of correlation functions. The first class. 



a = i^^jfj), b = {itvP = {VUP, and c = {^ivf^) 
involves operators which appear in the Hubbard basis 
and hence enter directly as self-consistent parameters in 
the closure for the Green function. In the second class, 
correlations of the form. 



(11) 



which corresponds to the bandwidth renormalization and 
g = Re(n^n^), arise between entities outside the Hub- 
bard basis. Both of these correlation functions must be 
decoupled and expressed in terms of correlation functions 
of the type a — c. Expressions appear in our earlier paper 
Q for the real part of the correlation function y. The 
important difference here is that now p is complex. Its 
imaginary part is given by 



where 



Imp = pi = 



(G22 + G12) — 



1 



1 - {n^o) 



(12) 



(Gn + G21) (13) 



and 



<5i 



[Rc(a + fe)Im(c + h) 



(n_,)(2-(n_,)) 
+Im(a + 6)Re(c + h)\ (14) 

where G12 = {iia^\a) = G21, Gn = (^fffL), and G22 = 

Noting that any correlation function of the form, 
{'4'm{'i)'4'}iij)) related to the Green function through 

n 



(2^) 



d2Mtje*'^-("---'-^)(l-/(w)) 



^(a,)„,„(k)<5(a;-£;,(k)) (15) 
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the complete set of equations for the correlation func- 
tions, a, 6, c, p, and g together with the equation for the 
Green function, Eq. (||), represent a closed set of inte- 
gral equations that must be solved self-consistently. In 
the absence of nearest-neighbour interactions, the terms 
containing a, b, and c are summed over nearest neigh- 
bour sites. Hence, the imaginary part cancels as a result 
of the alternating signs. Such is not the case for p, how- 
ever. Consequently, staggered current vorticity is stable 
if the closed set of equations admits a non-zero value of 
Imp = pi, the key signature of the broken symmetry. 
Shown in Fig. is the self-consistent solution for the 
lower Hubbard bands for U = 8t, n = 0.79, T = O.Olt 
and V = 0. 
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FIG. 1. Momentum dependence of the lower Hubbard band 
for (7 = 8t, n = 0.79, and T = O.Olt. From the spectral 
function shown in the inset, it clear that the spectral weight 
in the lower Hubbard band resides entirely on the curve with 
a fiat dispersion near (tt, tt). The spectral weight in the upper 
Hubbard band is also shown. 

The two bands in Fig. (|^) correspond to i?i(k) and 
i?2(k) — i?i(k -I- Q). As the spectral function at- 
tests (see inset Fig. (I)), only one of the bands has a 
non-zero weight as is expected in the absence of stag- 
gered currents. The spectral function shown in the inset 
-l/Trlm^^mn (G^^^ + ^^^)mn actuality consists of a 
series of 5-functions. Each peak has been given an artifi- 
cial width of 0.2t. 

To investigate the stability of staggered currents, we 
set pi to a fixed value, pin, and determine the self- 
consistent values of all correlation functions upon suc- 
cessive iteration. We then obtained the output value of 
Pi — Pout by using Eq. (|2|) . Staggered current vorticity 
is stable if p-^ = Pont or equivalently, A(pi„) = Pont/Pin = 
1 and dX/dpin < 0. As is evident from Fig. (||), A < 1 and 
it is a decreasing function of pin, proving that staggered 



currents are absent for U = 8t, V = 0.0 and T = O.Olt. 
This can also be seen from the inset of Fig. which 
contains the average energy per particle {{H)) as a func- 
tion of pi . The minimum value of the energy corresponds 
to pi = 0, further affirming the absence of staggered cur- 
rents. 
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FIG. 2. Staggered flux stability parameter A vs the input 
value (pin) for the imaginary part of the correlation function, 
p. The flux phase is stable if A = 1 and dA/p. < 0. Because 
A < 1, staggered currents are absent. The inset demonstrates 
that the energy per particle is a minimum when p is real, 
thereby implying an absence of the staggered current phase. 
Here U = 8t a.nd V = 0. 

Although the flux phase is absent for the parameters 
above, we can examine the tendency towards such order- 
ing as a function of density for a sufficiently small value of 
Pi such that A « A(0), the maximum value of A as shown 
in Fig. (^). Shown in Fig. (^) is the density depen- 
dence of the tendency towards staggered fiux ordering 
for pi = 0.001 for three values of the on-site Coulomb 
repulsion, U = 8t, 4t, 2t, V = 0.0 and T = O.Oli. 




filling (n) 

FIG. 3. Density dependence of the staggered current sta- 
bility parameter. Half-filling corresponds to n = 1. 
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As is evident, stable staggered currents are most likely 
at small values of U and as the density approaches half- 
filling as has been reported earlier in mean-field studies 
of the t — J model ||] . However, for the parameters of the 
Hubbard model that are relevant for the copper oxides, 
[/ » such ordering is absent at this level of theory. 
Going beyond the static approximation involves includ- 
ing quantum fluctuations. Such processes can only sup- 
press ordering. Hence, we conclude that for U ^ t, the 
Hubbard model (with V = 0) cannot sustain staggered 
flux ordering. 




-0.5 -0.25 0.25 
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FIG. 4. Staggered current stability parameter as a function 
of the nearest-neighbour Coulomb interaction, V for n = 0.79, 
U = 8t and T = Q.Olt. For V < 0, the stability parame- 
ter exceeds unity leading to an onset of staggered flux order. 
From the inset, the onset temperature when V — — 0.25t is 
T = 0.33t ^ mOK. 



Fig. (^)) of the staggered current is not entirely con- 
sistent with competing order in the underdoped regime. 
Hence, either staggered flux order is irrelevant to the cop- 
per oxides or at least nearest neighbour attractions are 
essential to their stabilization. An obvious source of 
such an attraction are local phonon modes as has been 
discussed recently Our work suggests (assuming the 
flux phase is relevant) that the strange phenomenology 
of the cuprates arises from cooperation between phonons 
and strong on-site Coulomb repulsions. 
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FIG. 5. Lower Hubbard energy bands and spectral function 
for n = 0.79, V = —0.25t and U = 8t. The presence of a gap 
along the (tt, 0) — (0, n) direction with a node at {n/2, n/2) is 
indicative of d^2_y2 order. 



What about nearest-neighbour Coulomb interactions? 
As is evident from Fig. (jij), the flux phase is desta- 
bilised even further for > 0. However, for < 0, 
A crosses unity and staggered flux ordering obtains. This 
is the primary result of this study. The inset shows the 
temperature dependence of A for = 0.0 as well as for 
V = — 0.25i. The onset temperature for the flux phase 
is T « 0.033i « 160K. For V = 0.0, A is relatively flat 
indicating further that even at T = 0, such ordering is 
absent. To examine the impact of the staggered current 
on the band structure, we recalculated the lower Hubbard 
band for V = -0.25t, n = 0.79, C/ = 8i and T = O.OIi. 
From Fig. (^), we see that a non-zero staggered current 
vorticity has opened a gap near the Fermi energy along 
the (tt, 0) — (0, tt) direction. The inset shows more clearly 
that the Fermi level lies in the gap near the (tt, 0) and 
(0, tt) regions. The position of the Fermi energy relative 
to the gap edge is doping dependent. 

Our work clearly shows that for the parameters rel- 
evant to the copper oxides, staggered current vorticity 
cannot be stabilized unless attractive interactions are 
present. Further, the apparent doping dependence (see 
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